#include"ppSpline.h"
#include<cmath>

// 计算两点间的距离
template<typename T>  
double distance_2pts(const vector<T>& pt1, const vector<T>& pt2)
{  
    if(pt1.size() != pt2.size()){
        return -1;  // 维数不匹配
    }  
    double sum = 0;  
    int n = pt1.size();
    for(int i = 0; i < n; i++){  
        sum += (pt1[i] - pt2[i]) * (pt1[i] - pt2[i]);  
    }  
    return sqrt(sum);
}

// 得到二维vector的第j列数
template<typename T> 
vector<T> get_column(const vector<vector<T>>& matrix, int j)
{  
    vector<T> column;  
    for (const auto& row : matrix){
        int N = row.size();  
        if(j >= 0 && j < N){  
            column.push_back(row[j]);  
        }
        else{  
            return {};  // 给定的列数超出范围
        }  
    }  
    return column;  
}  

ppSpline<3> solver(const vector<double> &knots, const Function &fun,
                    const boundary_condition boundary)
{
    int N = knots.size();
    vector<double> vals(N);
    for(int i = 0; i < N; i++){
        vals[i] = fun.f(knots[i]);
    }
    if(boundary == complete){
        return solver(knots, vals, boundary, 
                        fun.df(knots[0]), fun.df(knots[N-1]));
    }
    else if(boundary == natural){
        return solver(knots, vals, boundary);
    }
    else if(boundary == second){
        return solver(knots, vals, boundary, 
                        fun.d2f(knots[0]), fun.d2f(knots[N-1]));
    }
    else{
        cerr << "Error! 传入参数有误,请重试!" << endl;
        exit(-1);
    }
}

ppSpline<3> solver( const vector<double> &knots,
                    const vector<double> &vals,
                    const boundary_condition boundary)
{
    if(knots.size() != vals.size())
    {
        cerr << "Error!维数不匹配,请重试!" << endl;
        exit(-1);
    }

    int N = knots.size();
    vector<vector<double>> A;   //系数矩阵A
    vector<double> b;   //线性方程组右端项
    vector<double> u;   //A的上次对角线
    vector<double> l;   //A的下次对角线
    vector<double> d;   //A的主对角线  

    if(boundary == natural){
        A = get_diff_table(knots, vals);
        vector<double> b0 = get_diff_table_col(A, 3);
        for(int i = 0; i < N-2 ; i++) b0[i] *= 6;

        b.push_back(0.0);
        b.insert(b.end(), b0.begin(), b0.end());
        b.push_back(0.0);

        u.push_back(0.0);
        for(int i = 1; i <= N-2; i++){
            double lambda = (knots[i+1]-knots[i]) 
                            / (knots[i+1]-knots[i-1]);
            u.push_back(lambda);
        }

        for(int i = 1; i <= N-2; i++){
            double mu = (knots[i]-knots[i-1]) 
                        / (knots[i+1]-knots[i-1]);
            l.push_back(mu);
        }
        l.push_back(0.0);

        d.push_back(1.0);
        for(int i = 1; i <= N-2; i++){
            d.push_back(2.0);
        }
        d.push_back(1.0);       
    }
    else{
        cerr << "Error! 参数不匹配,请重试!\n";
        exit(-1);
    }
    
    vector<double> x = solu(l,u,d,b);
    vector<polynomial<3>> polys;
    for(int i = 1; i < N; i++){
        polynomial<3> S = cubic_0_2_val(knots[i-1], knots[i],
                                        vals[i-1], vals[i],
                                        x[i-1], x[i]);
        polys.push_back(S);
    }
    ppSpline<3> ppspline3(knots, polys);
    return ppspline3;
}

ppSpline<3> solver( const vector<double> &knots,
                    const vector<double> &vals,
                    const boundary_condition boundary,
                    const double &left_bound, const double &right_bound )
{
    if(knots.size() != vals.size())
    {
        cerr << "Error!维数不匹配,请重试!" << endl;
        exit(-1);
    }

    int N = knots.size();
    vector<vector<double>> A;   //系数矩阵A
    vector<double> b;   //线性方程组右端项
    vector<double> u;   //A的上次对角线
    vector<double> l;   //A的下次对角线
    vector<double> d;   //A的主对角线
    
    if(boundary == complete){
        vector<double> inter_pts;   
        vector<double> vals_of_pts;
        for(int i = 0; i < N; i++){
            inter_pts.push_back(knots[i]);
            vals_of_pts.push_back(vals[i]);
            if(i == 0){
                inter_pts.push_back(knots[i]);
                vals_of_pts.push_back(left_bound);
            }           
            if(i == N-1){
                inter_pts.push_back(knots[i]);
                vals_of_pts.push_back(right_bound);
            }    
        }

        A = get_diff_table(inter_pts, vals_of_pts);
        
        b = get_diff_table_col(A, 3);
        for(int i = 0; i < N; i++) b[i] *= 6;

        u.push_back(1.0);
        for(int i = 1; i <= N-2; i++){
            double lambda = (knots[i+1]-knots[i]) 
                            / (knots[i+1]-knots[i-1]);
            u.push_back(lambda);
        }

        for(int i = 1; i <= N-2; i++){
            double mu = (knots[i]-knots[i-1]) 
                        / (knots[i+1]-knots[i-1]);
            l.push_back(mu);
        }
        l.push_back(1.0);

        for(int i = 1; i <= N; i++){
            d.push_back(2.0);
        }
    }
    else if(boundary == second){
        A = get_diff_table(knots, vals);      
        vector<double> b0 = get_diff_table_col(A, 3);
        for(int i = 0; i < N-2 ; i++) b0[i] *= 6;
        b.push_back(left_bound);
        b.insert(b.end(), b0.begin(), b0.end());
        b.push_back(right_bound);

        u.push_back(0.0);
        for(int i = 1; i <= N-2; i++){
            double lambda = (knots[i+1]-knots[i]) 
                            / (knots[i+1]-knots[i-1]);
            u.push_back(lambda);
        }

        for(int i = 1; i <= N-2; i++){
            double mu = (knots[i]-knots[i-1]) 
                        / (knots[i+1]-knots[i-1]);
            l.push_back(mu);
        }
        l.push_back(0.0);

        d.push_back(1.0);
        for(int i = 1; i <= N-2; i++){
            d.push_back(2.0);
        }
        d.push_back(1.0);
    }
    else{
        cerr << "Error! 参数不匹配,请重试!\n";
        exit(-1);
    }
    

    vector<double> x = solu(l,u,d,b);
    vector<polynomial<3>> polys;
    for(int i = 1; i < N; i++){
        polynomial<3> S = cubic_0_2_val(knots[i-1], knots[i],
                                        vals[i-1], vals[i],
                                        x[i-1], x[i]);
        polys.push_back(S);
    }
    ppSpline<3> ppspline3(knots, polys);
    return ppspline3;
}

vector<ppSpline<3>> solver( const vector<vector<double>> &points,
                            const boundary_condition boundary )
{
    int N = points.size();   // 插值点个数
    int n = points[0].size();  // 每个点的维度

    vector<double> knots(N);
    knots[0] = 0;
    for(int i = 1; i < N; i++){
        double d = distance_2pts(points[i-1], points[i]);
        knots[i] = d + knots[i-1];
    }

    vector<ppSpline<3>> n_dim_ppSpline;
 
    if(boundary == natural){
        for(int j = 0; j < n; j++){
            vector<double> j_coord = get_column(points, j);
            n_dim_ppSpline.push_back(solver(knots, j_coord, boundary));
        }
    }
    else{
        cerr << "Error! 参数不匹配,请重试!\n";
        exit(-1);
    }

    return n_dim_ppSpline;
}

vector<ppSpline<3>> solver( const vector<vector<double>> &points, 
                            const boundary_condition boundary,
                            const vector<double> &left_bound, 
                            const vector<double> &right_bound )
{
    if(points[0].size() != left_bound.size()
    || points[0].size() != right_bound.size())
    {
        cerr << "Error!维数不匹配,请重试!" << endl;
        exit(-1);
    }

    int N = points.size();   // 插值点个数
    int n = points[0].size();  // 每个点的维度

    vector<double> knots(N);
    knots[0] = 0;
    for(int i = 1; i < N; i++){
        double d = distance_2pts(points[i-1], points[i]);
        knots[i] = d + knots[i-1];
    }

    vector<ppSpline<3>> n_dim_ppSpline;
 
    if(boundary == complete || boundary == second){
        for(int j = 0; j < n; j++){
            vector<double> j_coord = get_column(points, j);
            n_dim_ppSpline.push_back(solver(knots, j_coord, boundary,
                                        left_bound[j], right_bound[j]));
        }
    }
    else{
        cerr << "Error! 参数不匹配,请重试!\n";
        exit(-1);
    }

    return n_dim_ppSpline;
}

ppSpline<1> solver(const vector<double> &knots, const Function &fun)
{
    int N = knots.size();
    vector<double> vals(N);
    for(int i = 0; i < N; i++){
        vals[i] = fun.f(knots[i]);
    }
    return solver(knots, vals);
}

ppSpline<1> solver(const vector<double> &knots, const vector<double> &vals)
{
    if(knots.size() != vals.size())
    {
        cerr << "Error!维数不匹配,请重试!" << endl;
        exit(-1);
    }

    int N = knots.size();
    vector<polynomial<1>> polys;
    for(int i = 1; i < N; i++){
        polynomial<1> S = linear_0_val(knots[i-1], knots[i],
                                        vals[i-1], vals[i]);
        polys.push_back(S);
    }
    ppSpline<1> ppspline1(knots, polys);
    return ppspline1;
}

vector<ppSpline<1>> solver(const vector<vector<double>> &points)
{
    int N = points.size();   // 插值点个数
    int n = points[0].size();  // 每个点的维度

    vector<double> knots(N);
    knots[0] = 0;
    for(int i = 1; i < N; i++){
        double d = distance_2pts(points[i-1], points[i]);
        knots[i] = d + knots[i-1];
    }

    vector<ppSpline<1>> n_dim_ppSpline;
    for(int j = 0; j < n; j++){
        vector<double> j_coord = get_column(points, j);
        n_dim_ppSpline.push_back(solver(knots, j_coord));
    }
 
    return n_dim_ppSpline;
}
